EEG FREQUENCY ANALYSIS ON THE 
PDP LAB 8/E COMPUTER SYSTEM 



Lawrence Morrison Gorham 



DUDLEY KNOX LIBR'^R* 

NAVAL POSTGRADUATE SChOOi 
MONTEREY# CALIFORNIA 93940 




NAVAL POSTGRAE 

Monterey, California 



SCHOOL 




TH E3IS 

EEG FREQUENCY ANALYSIS ON THE ' 

PDP LAB 8/E COMPUTER SYSTEM 

by 

Lawrence Morrison Gorham 
September 1974 

•Thesis Advisor: G. K. Poock 



f 



Approved for public release; distribution unlimited. 



UNCLASSIFIED 



SECURITY CLASSIFICATION OF THIS PAGE (When Dmtm Entered) 



REPORT DOCUMENTATION PAGE 


READ INSTRUCTIONS 
BEFORE COMPLETING FORM 


1. REPORT NUMBER 


2. GOVT ACCESSION NO. 


3. RECIPIENT'S CATALOG NUMBER 


4. TITLE (and Subtitle) 

EEG Frequency Analysis on the PDP Lab 
8/E Computer System 


5. TYPE OF REPORT ft PERIOD COVERED 

Master's Thesis; 
September 1974 


6. PERFORMING ORG. REPORT NUMBER 


7. AUTHORf*) 

Lawrence Morrison Gorham 


8. CONTRACT OR GRANT NUMBER(e) 


9. PERFORMING ORGANIZATION NAME AND ADDRESS 

Naval Postgraduate School 
Monterey, California 93940 


10. PROGRAM ELEMENT. PROJECT, TASK 
AREA ft WORK UNIT NUMBERS 


11. CONTROLLING OFFICE NAME AND ADDRESS 

Naval Postgraduate School 
Monterey, California 93940 


12. REPORT DATE 

September 1974 


13. NUMBER OF PAGES 

42 


M MONITORING AGENCY NAME ft ADDRESS^// dltterent trom Controlling Office) 

Naval Postgraduate School 
Monterey, California 93940 


IS. SECURITY CLASS, (ol thlt report) 

Unclassified 


15*. DECLASSIFICATION/ DOWNGRADING 

schedule 


16. DISTRIBUTION STATEMENT (of thle Report) 

Approved for public release; distribution unlimited. 


17. DISTRIBUTION STATEMENT (of the abetract entered In Block 20, If different from Report) 


18. SUPPLEMENTARY NOTES 



19. KEY WORDS (Continue on reveree mlde It neceeemry end Identity by block number) 



EEG 

Fast Fourier Transformation (FFT) 
PDP Lab 8/E Computer System 



20. ABSTRACT (Continue on reveree elde If neceeemry mnd Identify by block number) 

This thesis describes the analysis, method and computer 
programs used to obtain the Fast Fourier Transformation 
(FFT) of an electroencephalogram (EEG) using a small labo- 
ratory computer like the PDP Lab 8/E. The EEG power spectrum 
was then computed from this transformation. The information 
contained in this thesis is intended to enable the user to 



DD 1473 

(Page 1) 



EDITION OF 1 NOV 65 IS OBSOLETE 
S/N 0 102-0 14- 6601 | 



UNCLASSIFIED 

SECURITY CLASSIFICATION OF 



THIS PAGE Detm Entered) 



1 



UNCLASSIFIED 

CuLIjmTY CLASSIFICATION OF THIS PAGEfHTi.n Data Ent'rtd) 



Block #20 Continued 



compute the Fourier coefficients of a set of data points 
or compute the power spectrum of a real waveform such as 
the EEG . 



DD Form 1473 (BACK) 
, 1 Jan 73 

S/N 0102-014-GG01 



UNCL AS SIFIED . „ 

SECURITY CLASSIFICATION OF THIS PAGEfWIion D or* Fnfrtd) 



2 



EEG Frequency Analysis on the 
PDP Lab 8/E Computer System 

by 

Lawrence Morrison Gorham 
Lieutenant, United States Navy 
B.S., Auburn University, 1967 



Submitted in partial fulfillment of the 
requirements for the degree of 



MASTER OF SCIENCE IN OPERATIONS RESEARCH 

from the 

NAVAL POSTGRADUATE SCHOOL 
September 1974 





!=: • ■?. - 
c- / 








DUDLEY KNOX LlBR Ry 
NAVAL POSTGRACUA I £ S ,n 
MONTEREY, CAL.EORN.A 



ABSTRACT 

This thesis describes the analysis, method and computer 
programs used to obtain the Fast Fourier Transformation (FFT) 
of an electroencephalogram (EEG) using a small laboratory 
computer like the, PDP Lab 8/E. The EEG power spectrum was 
then computed from this transformation. The information 
contained in this thesis is intended to enable the user to 
compute the Fourier coefficients of a set of data points or 
compute the power spectrum of a real waveform such as the 
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I. INTRODUCTION 



As cerebral cortex cells discharge electrically, a small 
potential is produced; when this potential is averaged with 
the other discharging cells it can be sampled with electrodes 
applied to the scalp. This weak voltage is amplified 100,000 
times by a low noise, high gain amplifier. For simplicity 
and economy, a paper and pen writing oscillograph is used to 
record the changing potential. This record of cortical po- 
tentials is called an electroencephalograph (EEG) . 

These waves are very complicated and analyzing them has 
been a problem since their discovery. Due to the wave-like 
nature, EEG frequency or spectrum analysis has been the 
subject of much research over the years. Large, high speed 
computers are being used to aid this work today. The pur- 
pose of this study was to develop fast Fourier techniques 
for spectrum analysis of EEGs using a small laboratory com- 
puter like the PDP Lab 8/E computer system. With the com- 
puter located in the laboratory, near subject and equipment, 
it was hoped that some kind of real time frequency analysis 
could be accomplished. 
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II. BACKGROUND 



A. PHYSIOLOGY 

The small potential measured in the EEG originates in 
the outer .1 inch of the brain called the grey matter. The 
working cells of the grey matter are the neurons, numbering 
perhaps to ten billion in all. Very little is known about 
how these cells operate to provide the many functions that 
the brain accomplishes. One well known and supported theory 
by John C. Eccles explains how the EEG waveforms are pro- 
duced [Ref. 10]. Each neuron is assumed to be an independent 
functioning unit, receiving stimuli from outside sensory 
devices and other neurons through dendrite fibers. When a 
particular number of stimuli are received, the neuron fires 
a pulse of electricity to other neurons via axon fibers. 
Receiving reverse polarity signals will cancel the effect of 
incoming stimuli and prevent the neuron from firing. In 
this way each neuron may be connected to as many as 100 
other neurons. 

Usually many thousands of neurons are stimulated simul- 
taneously. As they reach their excitation points and fire, 
other thousands of neurons receive these new stimuli. Thus 
the collection of stimuli moves through the brain step by 
step. Each stage in the path may be only a thousandth of 
an inch or it may be three or four inches in another part 
of the brain. After firing, each neuron has a deactivated 
period of several seconds, in which it cannot be stimulated 
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to respond. This delay, or refractory period as it is called, 
produces a wave-like quality to the stimuli collection as it 
moves along. Eccles believes that this explains the wave- 
like nature of cortical potentials as they filter out through 
the skull and scalp and are recorded on the EEG. 

Another important aspect of these neural pathways is that 
they frequently double back on each other. This causes a 
circular path to form, which provides a self -exciting fea- 
ture to the brain. These closed loops, with fixed distance 
and refractory periods , may account for the predominance of 
certain frequency waves. The prospect of neurons having 
feedback loops opens up a whole new level of complexity to 
the brain. In this way each neuron could be viewed as a 
very flexible amplifier, each one operating differently, 
depending on the kind of feed-back it is receiving. 

B. HISTORY 

A short history of the cortical potential must begin 
with Richard Canton, a British physiologist, who reported 
in 1875, "Feeble currents of varying direction pass through 
the multiplier when electrodes are placed on two points of 
the external surface (of the brain) , or one electrode on the 
grey matter, and one on the surface of the skull." At the 
time, Canton was examining cat, monkey and rabbit brains 
with a sensitive galvanometer using optical magnification. 

It is remarkable that he could distinguish such small po- 
tentials without the use of electronic amplification, which 



8 



wasn't developed for another fifty years. In 1876 the 
Russian, Danilevskey, discovered a change in the average 
potential of the cortex in response to acoustic stimuli. 

Only recently has this interesting aspect of brain activity 
been reinvestigated. In 1914 Cybulski recorded an apileptic 
seizure in a dog. The early attempts at recording cortical 
potential used photography, however, this was expensive and 
tedious. Pen recorders were not available for this use 
until 1940 [Ref. 8] . 

Recording cortical potentials in man lagged behind their 
demonstration in animals for many years. In man, both elec- 
trodes were attached to the scalp, causing the potentials 
to be greatly attenuated through the skull and scalp. The 
first recordings from the human brain were made in 1924 by 
Hans Berger in Germany. He worked and published much, ex- 
ploring many aspects of the use of the electroencephalogram. 
Some of these included physiological, neurological, psychia- 
tric and psychological applications-. He confirmed Canton's 
findings in animals and established that these waveforms 
change with age, with sensory stimuli and with chemical 
changes in the body. He was the first to record a major 
epileptic seizure in man. Due to the earlier work by 
Kaufman and Cybulski, it became known that epilepsy in 
animals caused abnormal waveforms. 

After Berger's first publication in 1929, many labora- 
tories took up this, work and knowledge in the field grew 
rapidly. After substantiating Berger's claims, the various 
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brain centers began to specialize in different research 
areas: Gibbs and Lennox in epilepsy, the Davises in normal 

subjects, Adrian, Hoagland and Jasper in the physiological 
mechanisms and Walter in brain lesions. The growth of 
electronencephalography has really paralleled the develop- 
ment of the machinery for recording the waves. In the 1930 's 
the electronic amplifier gave a tremendous boost to the 
field. In 1948 the transistor allowed the development of 
highly portable equipment. Modern instruments allow 16 
channels to be monitored simultaneously, each of which has 
linear response and accurately reproduces frequencies from 
1 to 500 Hz. A more complete history of electroencephalo- 
graphy may be found in Braizier [Ref. 4] . 

One of the many things that Hans Berger discovered about 
the human electroencephalograph was that the most normal 
pattern in adults, with eyes closed, was a frequency between 
8 and 13 Hz which he called alpha rhythm. These and other 
rhythms seem to originate in the cortex, however, they cease 
if the cortex is surgically disconnected from the thalamus. 
Alpha rhythm is most strong when the subject is relaxed 
with eyes closed. Opening the eyes in normal light signifi- 
cantly reduces alpha rhythm and similarly, opening the eyes 
in total darkness causes a temporary reduction. There is 
some evidence [Ref. 18] to support a positive correlation 
between high amplitude alpha rhythms and a passive or in- 
troverted personality structure as revealed by psychoanal- 
ysis. Similarly, low amplitude alpha rhythms were found 
prevalent in extroverts [Ref. 18] . 
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Delta activity is between 1 and 3.5 Hz and is generally 
associated with sleep. In 1954 Bremer achieved a strong 
delta state by disconnecting all the sensory pathways. This 
is logical since sleep is most easily induced when sensory 
stimulation is reduced to a minimum. It must be emphasized 
however, that even though no light, sound, etc. is sensed, 
one does not sleep unless he is "tired," i.e., delta state 
and sleep are not the same. Theta activity, in the 4 to 7 
Hz range, is usually associated with normal children [Ref. 4]. 
The percentage of theta decreases significantly with age in 
normal subjects. Beta activity, with a frequency of 14 to 
30 Hz, is generally associated with increased anxiety levels. 
The incidence of beta is increased in subjects under-going 
alcohol withdrawal. Beta activity can easily be blocked in 
the motor center associated with a particular movement [Ref . 

4] . 

The four basic waves, alpha, delta, theta and beta, span 
the frequency spectrum of cortical potentials. Another pat- 
tern of interest is the mu rhythm, which is characterized 
by a sharp peak in negative half -cycle, rounded positive 
half -cycle and 7 to 11 Hz. Mu rhythm is reduced not only 
by overt movement of a limb but also in some cases when 
the subject is contemplating movement [Ref. 4, 8]. 

Except for beta waves, all the waveforms previously 
described tend to be reduced by some kind of stimulation. 

In contrast, lambda waves, a saw-toothed waveform, are best 
seen when the subject is actively looking at something of 



11 



interest to him. This waveform usually occurs when the 
subject is scanning a complex pattern or picture. Paradox- 
ically, random waves of the same polarity and shape occur 
in some subjects during light sleep. A thorough analysis 
of the various waveforms in psychologically normal and 
abnormal subjects may be seen in Thorp and Heyman [Ref. 24] 
and Kiloh and Osselton [Ref. 21] . 

Interpretation of the EEG by visual inspection of the 
record depends in a large part of recognizing common pat- 
terns and their association with psychologically normal or 
pathological conditions. Unfortunately these visual methods 
are unable to differentiate between clinical conditions 
having similar EEG's. Several different methods for anal- 
yzing the EEG are being tested. One of the first attempts 
was to get a quantitative grasp on the amplitude of the 
waves. Walter [Ref. 4] joined the peaks and troughs with 
two lines and then measured the average distance between 
the lines. With the advent of amplifiers, electronic inte- 
gration solved this problem. 

In looking for better ways to specify cortical poten- 
tials, many forms of frequency analysis have been tried. 

As early as 1932, Dietsch suggested that these potentials 
might be composed of sine waves and he outlined the procedure 
whereby Fourier analysis could be used to determine the fre- 
quency components present [Ref. 9] . 

The first frequency analysis of an EEG was accomplished 
at the Cincinnati College of Medicine in 1944, where the fre- 
quencies were determined by manually counting the number of 
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complete waves in each second [Ref. 12]. In the late for- 
ties and fifties, much work was stimulated by Gibbs and 
Grass [Ref. 13] when wide-band electronic filters were used 
to separate the different frequencies [Ref. 19, 22]. Filter 
devices are commonly used in this way, even today [Ref. 25]. 
An excellent review of all work on EEG analysis techniques 
up through the fifties may be found in Reference 5. 

Spectrum analysis of physiological data using the dis- 
crete Fourier transform was accomplished by Halberg in 1961 
[Ref. 14, 23]. Walter was first to use the discrete trans- 
form on cortical potentials in 1963 [Ref. 26, 28], Most of 
the quantitative frequency work being done on EEG's today 
is with the aid of large computers in which the discrete 
Fourier transform provides the conversion [Ref. 11, 15, 17]. 
The U. S. Navy is also entering this research area, at the 
Navy Medical Neuropsychiatric Research Unit in San Diego, 
California [Ref. 15] . In 1973 the Fast Fourier Transform, 
a faster version of the discrete Fourier transform, came 
into use on the EEG [Ref. 2] . 
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III. FOURIER ANALYSIS 



A. FOURIER TRANSFORM 
1 . Continuous 

An introduction to Fourier analysis must begin with 
the general form of the Fourier series. It is 



f(t) = A^ + £ , A„ cosnoot + „£., B sinnoot (1) 

v ' on=ln n=ln v J 



where go = 2ir/T and T is the period. Any waveform can be 
represented by determining the coefficients A , A^ and 
and summing all the terms of this series. The coefficients 
give the relative contribution of each frequency that is 
present in the signal. These are harmonics or integer 
multiples of the fundamental frequency 1/T. Fourier de- 
veloped a set of equations from (1) and solved them simul- 
taneously. The result 
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provides a method for finding the harmonics when the function 
is known. 

2 . Discrete 

Frequently, however, the complete function is not 

known and only a discrete number N of samples over a finite 
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period T is known. In this situation, an approximation to 
the continuous function may be written 



„ , „ n n 

f(t) = A^ + 2 E, A cosmoot + 2E, B sinmwt (5) 
^ J o m=l m m=ln v J 

where n = N/2 for N even and n = N+l/2 for N odd. The ex- 
pressions for the coefficients now become 
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( 6 ) 

(7) 



These functions are easily computed from the set of 
sample data points; however, the time required is sometimes 
prohibitive. It will be noted that this procedure requires 
N 2 trigonometric operations, since each of N coefficients is 
computed from N terms of a series. These computations, in 
addition to numerous other arithmetic operations, require 
considerable computer time, especially on small laboratory 
computers like the PDP Lab 8/E. 



B. PROBLEM AREAS 

There are three problems commonly encountered when using 
the discrete Fourier transform. These are: aliasing, leak- 

age and the picket-fence effect. Aliasing occurs when a 
higher than expected frequency is sampled, but appears in 
the spectrum as a lower frequency. Suppose an unknown sig- 
nal is being sampled five times at a rate of 100 Hz. If 
an 80 Hz component is present, it will appear in the spectrum 
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as 20 Hz and there is no way of distinguishing this pseudo- 
frequency from the real one. Nyquist proposed the solution, 
in that, the sampling rate should be twice the highest fre- 
quency present [Ref. 16]. This can' be accomplished in an 
unknown signal by passing it through a lowpass filter set 
at half the sampling rate. A lowpass filter was used to 
prevent aliasing in this study. 

The problem of leakage occurs in all Fourier analysis 
using finite sampling records. The record, or the total 
length of the sample, can be viewed as a window. When this 
window is transformed, without the sample data present, the 
resulting spectrum contains not only the fundamental fre- 
quency but also spurious peaks called sidelobes. When the 
window with its sample is transformed, these sidelobes tend 
to leak into the spectrum. Because of the relatively small 
sidelobes, in the frequency range of interest, leakage was 
not a problem in this study. 

The picket-fence effect is a problem when a pure sinu- 
soid occurs between two adjacent harmonics. If it is ex- 
actly halfway between computed harmonics, the signal is 
reduced to .637 of the original in both harmonics. The 
worst situation occurs when there is a frequency beti^een 
each pair of adjacent harmonics. In this case an erroneous 
ripple, or picket-fence, appears in the spectrum. Pure 
sinusoids were assumed not to be present in EEG’s and thus 
the picket-fence effect was ignored. Bergland discusses 
these problems as well as the mathematics associated with 
them in Reference 3. 
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IV. FAST FOURIER ANALYSIS 



A. FAST FOURIER TRANSFORM 

In 1965 James Cooley and John Tukey published an impor- 
tant paper that described a much faster method to compute 
the discrete Fourier transformation, called the Fast Fourier 
Transform (FFT) [Ref. 7] . The procedure adopted for use on 
high speed computers has become known as the Cooley-Tukey 
Algorithm. A general outline of the proof of this algorithm 
follows . 

A time series Xj, containing N even samples, can be par- 
titioned into two functions Y r and Z r , each containing N/2 
points. Let Y r be composed of the even-numbered points 
(x q ,X£ , . . . ,x^, 2 ) an ^ the odd-numbered points (x^ ,x^ . . .x^ _^j ? ) . 
Equations (6) and (7) may be combined using complex number 
notation, so that 



n- 1 

A = 
m r=0 

where W = and 

N/2 points, each has 
coefficients are 



X W mr 
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j = Z 7 !- 

a discrete 



Since Y 
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Fourier 



( 8 ) 

and Z are a set of 
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transform whose 



N/2-1 



B = Z n 
m r=0 



Y W 1 
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mr 



( 9 ) 



C = 1 zi 1 Z W m ( 2r+1 ) 
m r=0 r 



N 



( 10 ) 



with m = 0,1,...,— -1. Combining equations (9) and (10), to 
get odd and even points, gives 
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N/2-1 
A m = R=0 




m N/2-1 

+ w m z z w mr 

W R= 0 r 




m 



( 11 ) 



A + N/2 = B - W m C 
m m 



m 



( 12 ) 



The details of this proof will be found in Reference 6. 

If X is an even power of two, then Y and Z also have 
r r r 

an even power of two points and each can be partitioned again 
and again until only one point remains in each set. The 
Fourier transform of a single point is just that point again; 
the entire transformation can be accomplished using equations 
(11) and (12) repeatedly at each partition. The requirement 
that the number of points N be a power of two assures that 
after log 2 N partitions there will be only one point in each 
set. In this way, the entire set of points can be translated 
without carrying out the manipulations of equations (6) and 
(7) . Since the FFT is derived from the discrete Fourier 
transform, the FFT maintains the same precision as before. 

B. FLOWPATH 

An example of the FFT for N = 8 will be illustrative. 
Figure 1 shows the flow diagram and the configuration after 
the first partition. Using Tukey's notation, each dot rep- 
resents a variable and the factors entering the dot are 
additive. When a multiplication is required before adding, 
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Figure 1. Fast Fourier Transform Flov; Graph, First Partition. 
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the quantity will be beside the factor to be multiplied. 

The set of for j = is computed from equations 

(11) and (12) . The quantity near the bottom right dot is 
Ay = + W 7 ^. Figure 2 shows all quantities after two 

partitions and Figure 3 shows the complete transformation 
broken down into only multiplication and addition operations. 

It may be noted in Figure 3 that the data points must be 
rearranged according to the bit reversal rule before the 
algorithm can be carried out. Bit reversal is accomplished 
by representing the data point position number in binary 
and reversing the bits about the center. Position (3) is 
written (0 1 1) and is reversed to (1 1 0) , which is (6) . 

Thus (0) stays a (0) , (1) becomes a (4) and so on. More 
information about FFT may be found in References 3 and 6. 

C. PROGRAM 

The computer program FAST was written to accomplish the 
FFT on the PDP Lab 8/E computer system. A complete listing 
of FAST may be found in Appendix A. The program is written 
in BASIC, one of the commonly used languages on the PDP Lab 
8/E. FAST receives data points as inputs and outputs the 
Fourier coefficients. 

In order to operate the FFT, read FAST into core using 
normal methods, and then type RUN. The program will type 

ENTER DATA POINTS 
? 

and go into a wait loop. Type in 64 data points in their 
proper, order. The number of data points may be changed to 
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Figure 3. Fast Fourier Transform, Flow Graph, Last Partition. 
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any power of two, but cannot exceed 64 due to core limita- 
tions, which is 8,000 words in the PDP Lab 8/E. FAST will 
translate the data points and type the Fourier coefficients. 
The time required to convert 64 points is approximately 25 
seconds. The power associated with each harmonic can be 
obtained by squaring and adding the coefficients at each 
frequency. 

D. EXAMPLE 

As an example to illustrate the use of FAST, the square 
wave shown in Figure 4 was translated and the results com- 
pared with the coefficients from the continuous Fourier 
transform; Table I shows the result of a FAST translation 
of the square wave along with the coefficients. To deter- 
mine the A n coefficients using the continuous transform, 
equations (2) and (3) are integrated from -T/2 to T/2. 

Since f(t) is an odd function, the A n terms are all zero. 
Using equation (4), however, and integrating over the same 
range, it is easily shown that 



B 



n 



2 



- 2CosnTr 
mr 



(13) 



To compare with FAST, the same 32 harmonics were computed 
using equation (13) (see Table II) . The answers from FAST 
are within two decimal places of the answers computed from 
the continuous Fourier transform. 
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Figure 4. Square Wave. 




Figure 5. EEG Wave and Timing Pulse. 
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Table I. FAST Coefficient 
of Figure 4. 
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B(N) 

-8.5070 59E + 37 
1.273241 
-4. 74319 IE -09 
.4244135 
-2 . 371595E- 09 
.2546481 
-1. 581064E-09 
. 1818915 
-1 . 185798E-09 
.1414712 
-9 . 486382E- 10 
. 1157491 
-7. 905318E-10 
. 09794159 
-6. 77598 7E-10 
.08488271 
-5.928989E-10 
.07489651 
-5 . 270212E- 10 
. 06701267 
-4 . 743191E- 10 
.06063051 
-4 . 311992E- 10 
. 05535829 

- 3 . 952659E- 10 

.05092963 

- 3 . 648608E- 10 

.04715706 
-3. 387993E-10 
. 04390485 
-3. 162127E-10 
.04107228 



Continuous Fourier Transform of B 
for Square Wave of Figure 4. 



Coefficients 
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V. SPECTRUM ANALYSIS 



A. SPECTRUM 

In order to do frequency analysis of EEG waveforms, a 
computer program called SPECTRUM was written which incor- 
porated FAST (Appendix A) . The problems encountered were 
many. First, a relatively noise-free signal had to be pre- 
sented to the PDP Lab 8/E analog- to-digital (A-D) converter. 
Then the signal had to be sampled at the right rate and 
duration. Finally, the data had to be delivered to FAST in 
a compatible format for translation. The power spectrum was 
chosen as the output format since it has more intuitive ap- 
peal than some other measure such as voltage or current. 
According to Parseval's theorem, the power at each frequency 
can be computed by squaring the two coefficients and adding 
them together. This was done and the power spectrum is dis- 
played on the teletype as percent of total power at each 
harmonic . 

B. APPARATUS 

The Beckman type RM Dynograph Recorder was used to pick- 
up and amplify the cortical potentials. Using sodium chlo- 
ride paste, silver electrodes were applied to the subject's 
scalp at positions and T^. A standard Beckman A-C 
coupler was used with the TIME CONSTANT switch set to .03 
and HIGH FREQ switch on 3. Depending on the signal strength, 
the preamplifier was set to either 1 mv/cm or .5 mv/cm and 
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the power amplifier set for a multiplication of .02. This 
adjustment was made in order to maintain a relatively con- 
stant signal level from the power amplifier. The FILTER 
switch was set on 3, which is a 60 Hz lowpass filter. 

The EEG was recorded on one track of a Hewlett-Packard, 
four track, instrumentation recorder on 1.0 mil magnetic 
tape at 3-3/4 in/s tape speed. A seven ms timing pulse was 
recorded on another track from an Ortec pulse generator. 

The timing pulse made it possible for the computer to find 
a particular place on the tape. The two tracks were then 
connected to the A-D converter of the PDP Lab 8/E computer. 

As an alternative hookup, the A-D converter can be connected 
directly to the Beckman and sample the EEG as the subject 
produces it. 

C. SAMPLING RATE 

With the proper signal on channel 1 of the A-D converter, 
SPECTRUM samples and translates it to the power spectrum. 
According to Nyquist's rule, the sampling rate must be twice 
the maximum frequency. The maximum frequency passed by the 
Beckman filter is 60 Hz, so the sampling rate must be at 
least 120 Hz - . In order to give a little extra margin, 128 Hz 
was chosen so that the harmonics would come out in even 
multiples of two. FAST can only convert 64 points and sam- 
pling at 128 Hz, the period (T) of the sample is % second. 

The reciprocal of the period (1/T) determines the increment 
of the harmonic. In this case 1/T =• 2, thus each harmonic 
is exactly 2 Hz apart. 
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D. EXAMPLE 



After loading SPECTRUM and typing RUN, the program sam- 
ples, converts the data, and types the power spectrum. 

Table III shows an example of one SPECTRUM run, taken from 
a subject who was relaxing with his eyes closed. The pro- 
gram starts sampling when it receives a numbered timing 
pulse on channel three. The EEG waveform from which this 
example was taken, along with its timing pulse, is shown 
in Figure 5. SPECTRUM uses the timing pulses to locate the 
desired sample on the tape. 

E. EXPERIMENT 

By choosing the proper segment with the timing marks and 
replaying the tape 20 times, SPECTRUM was used to scan a 
length of tape arbitrarily set for 10 seconds and produce 
the power spectrum for each half second record. The results 
of this experiment were reduced to the graph of Table IV. 

It plots frequency against time, with each number represent- 
ing the percentage of total power for a given frequency and 
time. If the table is viewed as a three dimensional graph, 
where percentage power is shown coming out of the page, it 
is easy to visualize a ridge running up the 10 Hz line. The 
subject was again resting with eyes closed, so the prepon- 
derance of alpha frequency (8, 10, 12 Hz) was understandable. 
The frequency also appeared to come in bursts. It was first 
on for about three seconds, then reduced for one second, then 
back on for a second, and so on for the length of the tape. 
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RUN 



? 



FREQUENCY 

HZ 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 



POWER 

% 

16.2738 

1.157298 

. 07388068 

4.684211 

5.177504 

3.729045 

49.26551 

5.911107 

2.511872 

3.593473 

. 7096686 

.6507585 

.9833005 

.9566361 

.6030759 

.538035 

.2399192 

.2305483 

. 3325844 

.2973473 

.2739161 

. 1222699 

.2039538 

. 178039 

. 1704812 

. 1763688 

. 1824439 

. 1651099 

. 1600951 

. 1526896 

. 1463384 

. 148733 



READY . 



Table III. 



A Typical SPECTRUM Run from a Half Second Sample 
of an EEG . 
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0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 

FREQUENCY (Hz) 



Table IV. Power Spectrum derived from an EEG Over a Period of 
Ten Seconds. Each Number Represents the Percentage 
Power at each Frequency for the Half Second Sampled. 
Zeros are omitted. 
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Sometimes this alpha burst effect can be detected right on 
the EEG. 

There were two places where strong delta frequency was 
present. Prepresented on the graph by the 2 Hz line, delta 
was present in the beginning and again around the five 
second point. This normally indicates a very light stage 
of sleep, where the subject is drifting in and out of con- 
sciousness. Under questioning the subject indicated that 
he could have been asleep several times during the recording 
session. The graph also shows the rapidity with which fre- 
quency changes [Ref. 2, 11, 15 and 27]. One example of this 
is from the first half second to the one second point, where 
the 10 Hz power changed from 3% to 39%. Rather than changing 
smoothly, it was as if one frequency stopped and another 
started . 

F. NOISE 

Either because of a noisy amplifier or possibly a dirty 
tape, it may become necessary at some time to analyze a 
noisy signal. This will become known when repeated spectra 
from the same point show a variance. In order to reduce 
this variance, a method of averaging was devised. Instead 
of sampling at 128 Hz, a new rate of 1280 Hz was used and 
the sample means were computed every 10 points. The sample 
period was kept at h second, so a total of 64 sample means 
were computed each run. This reduced the variance by a fac- 
tor of /To and providing the signal' was not too noisy, the 
results were satisfactory. Thus to obtain more information 
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from the noisy EEG, it may become necessary to sample more 
frequently than according to Nyquist’s law. 

Increasing the sampling rate is not without its short- 
comings. Firstly, it takes longer to process the added data. 
Changing to a 10 point average increases the computing time 
from 25 to 40 seconds. Secondly, part of the data must be 
stored while the averages are being computed. The new pro- 
gram requires about 200 more storage positions, which is at 
the maximum capability of the PDP Lab 8/E with 8,000 word 
memory. Most spectral analysis work being done on EEG's 
today use a 200 Hz sampling rate [Ref. 2, 15]. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



A real time analyzer for the EEG was one of the goals 
of this study. This has only partly been realized since 
SPECTRUM requires 25 seconds to convert a half second sam- 
ple. Since the EEG frequency pattern changes so rapidly, it 
is necessary to convert more than a half second; perhaps as 
much as 10 seconds in order to get a good representation of 
the changing nature of the wave. This means- that a much 
larger computer is needed to do real time spectral analysis 
of EEG ' s . 

As was shown using off-line techniques, SPECTRUM is a 
useful tool for spectral work. Using an instrumentation 
recorder in conjunction with some kind of pulsing device 
for timing, many different forms of frequency work may be 
accomplished such as autocorrelation, cross spectral analysis 
and digital filtering. SPECTRUM can also be used to compute 
power spectra from other kinds of continuous waveforms such 
as electrocardiograms and electromyograms. 

The author believes that this area of study is an impor- 
tant one and should be extended in the future; however, 
there are limitations. The major limitation to this system 
resides in the lack of adequate data storage, both in the 
small core size, which is about 4,000 words with BASIC in 
core, and in the lack of a random access device such as a 
disk pack. Of these two, the disk pack would be most useful 
for spectrum analysis work since large amounts of data from 
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the A-D converter could be stored rapidly and processed as 
required. The flexibility of the entire system would be 
considerably enhanced because the disk pack could also be 
used to break up large programs which are now too big for 
the computer and then process them in a segmented fashion. 
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APPENDIX A 



SPECTRUM 

1 N=64 

2 Nl=63 

5 DIM F(l,63) ,W(1 ,31) ,T(1) 

8 SET RATE 6,558 
10 DIM PI (300) 

15 REAL TIME PI, 1,0, 900 
18 IF ADC (3) <= . 9GO TO 18 
20 ACCEPT 
30 FOR 1=0 TO 63 
35 FOR J=1 TO 14 
40 F(0,1)=F(0,I) +ADB (0) 

45 NEXT J 
50 NEXT I 
58 PRINT 
100 FOR 1=0 TO 31 
110 READ F(1,I) 

115 F(1,63-I)=63-F(1,I) 

120 NEXT I 

158 FOR 1=0 TO N1 

162 IF F (1 , I) =1 THEN 169 

163 X=F (0 , 1 ) 

164 J=F (1,1) 

165 F (0 , I) =F (0 , J) 

166 F (0 , J) = X 

167 F(1 , I) =F (1 , J) 

168 F (1 , J) = J 

169 NEXT I 

194 FOR 1=0 TO N1 

195 F (1 , 1 ) =0 

196 NEXT I 

210 W( 0 , 0)=1 

220 W(0 , 1) =COS (6 . 2 83185/N) 

230 W(1,1)=-SIN(6. 283185/N) 

240 IF N<=4 THEN 295 
250 FOR 1=2 TO N/2-1 
260 11=1^1 

270 W(0,I)=W(0,1)*W(0.,I1) -W(l , 1) *W(1 , II) 
280 W(1,I)=W(1,1)*W(0,I1)+W(0,1)*W(1,I1) 
290 NEXT I 
310 J=1 

320 FOR M=1 TO N 
330 J=J*2 

340 IF J>=N THEN 370 
350 NEXT M 
370 FOR K= 1 TO M 
380 J2=2t (K- 1) 

390 1 2=N/ (J2* 2) 
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400 FOR 1=0 TO 12-1 

401 T(0)=F(0,Y) 

410 FOR J=0 TO J2-1 
500 X=N* I / 1 2 + J 
510 Y=X+ J2 
520 Z= J* I 2 

530 T(0)=W(0,Z)*F(0,Y) -W( 1 , Z) *F ( 1 , Y) 

540 T (1) =W(1 , Z) *F (0 , Y) +W(0 , Z) *F (1 , Y) 

550 F(0,Y)=F(0,X) -T(0) 

560 F(1,Y)=F(1,X)-T(1) 

570 F(0,X)=F(0,X)+T(0) 

580 F(1,X)=F(1,X)+T(1) 

590 NEXT J 
600 NEXT I 
610 NEXT K 
690 A=0 

700 FOR 1=0 TO N/2-1 

710 F(0,I)=(F(0,I)+2) + (F(l,I)+2) 

712 A=A+F (0,1) 

720 NEXT I 

728 PRINT "FREOUENCY POWER" 

730 FOR 1 = 0 TO 'n/2-1 
732 A2=F (0 , I ) /A 
750 PRINT 1*2, A2 
760 NEXT I 
800 FOR 1=0 TO N1 
810 F (0 , 1 ) = 0 
820 NEXT I 

900 DATA 0,32,16,48,8,40,24,56,4,36,20,52,12,44,28,60,2 

901 DATA 34,18,50,10,42,26,58,6,38,22,54,14,46,30,62,1 

998 RESTORE 

999 GO TO 1 



FAST 

1 N=64 

2 Nl= 6 3 

5 DIM F (1 , 63) ,W ( 1 , 31) ,T (1) 
30 FOR 1=0 TO N1 
35 INPUT F (0 , I ) 

49 F (1 , I ) =0 

50 NEXT I 
58 PRINT 

100 FOR 1=0 TO 63 
110 READ F ( 1 , 1 ) 

120 NEXT I 
130 FOR 1=0 TO N 1 
140 I 1= I *64/N 
150 F(1,I)=F(1,I1) 

155 NEXT I 

158 FOR 1=0 TO N1 

162 IF F (1 , 1 ) = I THEN 169 
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163 

164 

165 

166 

167 

168 

169 

194 

195 

196 

210 

220 

230 

240 

250 

260 

270 

280 

290 

310 

320 

330 

340 

350 

370 

380 

390 

400 

401 

410 

500 

510 

520 

530 

540 

550 

560 

570 

580 

590 

600 

610 

698 

699 

700 

710 

720 

730 

736 

738 

739 

740 

760 

761 



X=F(0,I) 

J=F(1,I) 

F (0 , I ) =F (0 , J) 

F(0 , J) =X 
F(1,I)=F(1,J) 

F(1,J)=J 
NEXT I 

FOR 1=0 TO N1 
F (1 ,I)=0 
NEXT I 
W ( 0 , 0)=1 

W(0 , l) = COS (6 . 283185/N) 

W(1,I)=-SIN(6. 283185/N) 

IF N<=4 THEN 295 
FOR 1=2 TO N/2-1 
11 = 1-1 

W(0,I)=W(0 ,1) *W(0,I1) -W(1,1)*W(1,I1) 
W(1,I)=W(1,1)*W(0,I1)+W(0,1)*W(1,I1) 
NEXT I 
J=1 

FOR M= 1 TO N 
J= J*2 

IF J>=N THEN 370 
NEXT M 
FOR K=1 TO M 
J2=2t (K- 1) 

I 2=N/ (J2*2) 

FOR 1=0 TO 12-1 
T(0) = F(0,Y) 

FOR J=0 TO J2-1 
X=N* I/I2+J 
Y=X+J 2 
Z=J* 1 2 

T(0)=W(0,Z)*F(0 , Y) -W(1,Z)*F(1,Y) 
T(1)=W(1,Z)*F(0,Y)+W(0,Z)*F(1 ,Y) 

F (0 , Y) =F (0 , X) -T (0) 

F(1,Y)=F(1,X)-T(1) 

F(0,X)=F(0,X)+T(0) 

F(1,X)=F(1,X)+T(1) 

NEXT J 
NEXT I 
NEXT K 
A=0 

PRINT "HARMONIC A(N) B (N) " 

FOR 1=0 TO N/2-1 
PRINT I,F(0,I)/32,F(l,I)/32 
NEXT I 

FOR 1=0 TO N/2-1 
Vl=2* I/N 
V2=F(0,I)/A 

IF V2<1. 000000E-03GO TO 760 
PRINT 1*2 ;V2 
NEXT I 
STOP 
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800 FOR 1=0 TO N1 
810 F (0 , I) =0 
820 NEXT I 

900 DATA 0,32,16,48,8,40,24,56,4,36,20,52,12,44,28,60,2 

901 DATA 34,18,50,10,42,26,58,6,38,22,54,14,46,30,62,1,33 

902 DATA 17,49,9,41,25,57,5,37,21,53,13,45,29,61,3,35,19 

903 DATA 51,11,43,27,59,7,39,23,55,15,47,31,63 

998 RESTORE 

999 GO TO 1 
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